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Abstract. We consider the evolution of a parton system which is formed at the 
central rapidity region just after an ultrarelativistic heavy ion collision. The evolution 
of the system, which is composed of gluons, quarks and antiquarks, is described 
by a relativistic Boltzmann equations with collision terms including radiation and 
retardation effects. The equations are solved by the test particle method using Monte- 
Carlo sampling. Our simulations do not show any evidence of kinetic equilibration, 
unless the cross sections are artificially increased to unrealistically large values. 
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1. Introduction 

The possible formation and evolution of a quark-gluon plasma (QGP) in relativistic 
heavy-ion collisions has been an important and very active research subject in recent 
years. Besides the theoretical studies of possible signatures of QGP formation, many 
investigations have been concerned with space-time dynamics of the QGP, its formation 
and final decay into hadrons. Perhaps the most common description of the evolution of 
the QGP has been based on relativistic hydrodynamics [[[], [|. This approach, however, 
does not address many important issues, in particular, the initial formation and local 
equilibration of the QGP. To explore this issue theoretically, a microscopic description 
of the parton system is required, which allows one to follow the phase-space evolution 
in detail. One such approach is the parton cascade model |3|, |], || ||. 

An important conceptual problem faced by such microscopic models is that they 
require a very detailed specification of the initial state. Due to our ignorance of 
nonperturbative processes in QCD, we do not know how to extract the required 
information from the parton wavefunctions of the colliding nuclei. Since the original 
formulation of the parton cascade model, important progress in our understanding of 
the parton structure of large nuclei at small Bjorken-x physics has been made 0. 
There is good reason to believe that the parton distribution at the beginning of a 
relativistic heavy ion collision can be described by means of semiclassical methods. 
This improvement of our understanding the initial state has motivated a resurgence 
of theoretical studies of the microscopic evolution of the matter formed in the central 
rapidity region in a nuclear collisions || H . 

To date, those studies have mostly focused on elastic collisions among partons, 
although the parton cascade model was originally formulated to include radiative 
processes as time-like and space-like branchings in the leading logarithmic approximation 
@. It has been pointed recently by Baier et al. || that the inelastic gg — > ggg process 
is the essential driving force toward the kinetic equilibration of the parton system. The 
importance of this process for the chemical equilibration of the QGP was already studied 



and pointed out by Biro et al. [10 



There is also a serious defect in many numerical implementations of the cascade 
model: the superluminal propagation of interactions [11], [5], |J . This becomes important 
whenever the hard-sphere picture of scattering interactions is applied without proper 
incorporation of the retardated nature of interactions in quantum field theory. 

In this article, we develop a cascade code which includes the process(gg — > ggg) 
as well as a relativistic picture of scattering that eliminates the superluminal problem. 
We assume a boost-invariant initial state and analyze the time evolution of the parton 
system. Because we do not include a model of the hadronization process, our description 
is useful only for a few fm/c after the onset of the reaction. 

In section 2 we present the equations of motion describing the evolution of a 
relativistic system composed of quarks, antiquarks, and gluons. In section 3 we explain 
the cascade code, especially the choice of the initial state and the parton dynamics. Our 
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results and conclusions are discussed in section 4. 



2. Equations of Motion 

We consider here a system, which consists of gluons, quarks (u, d, s) and their antiquarks, 
formed just after the impact of a heavy ion collision. The typical de Broglie wave 
length of a particle in the system, A = h/pc, is much shorter than the average distance 
between particles even at relatively high density so that we assume those particles can 
be considered as relativistic classical particles. We thus assume that the semi-classical 
but relativistic approach is applicable to this highly energetic and dense parton system. 
This may not be true for some particle whose energy is less than 1 GeV, but we limit 
our consideration to observables for which the particle picture is expected to be a good 
approximation. 

The system which we are studying can be described by the Lorentz-covariant 
Boltzmann equation jDJ, 

p u d ll f(x,p) = C(x,p,t) 1 (1) 

where f(x,p) is a phase space distribution function and the right-hand side denotes the 
collision terms. We assume that there is no external force acting on the system and the 
self-generated forces of the system are fully described by collisions. 

Gluons and quarks (u, d, s) and their antiparticles (u, d, s) are independent 
constituents so that the system is a mixture of 7 different kinds of partons. These 
particles can scatter on each other elastically and they can produce new particles. We 
will consider here only the dominant processes among many possible ones: gg «-> gg, 

gq <-> gq, gq <-> gq, gg -> 999, qq <-> qq, qq <-> qq, qq <-> qq, gg <-> qq- 

Including only those processes among partons, the Boltzmann equation for gluons 
takes the form: 

p"dMx,p)= / 2 / 3 /Jw„[/ 9 (3)/ ff (4)-/ g (l)/ ff (2)] 

+ / / / W gq ^ gq [f g (3)f q (A) - f 9 (l)f g (2)} 

+ 1 1 1 1 \w gg -> ggg [f g {±)f g {h) - / s (l)/ s (2)] 

+ / / /^™[A(3)/ 9 -(4)-/ ff (l)/ g -(2)] 

lLLl [Wq ^" mm - W 99^Ml)W)l (2) 

where we used the abbreviated notations f g (i) = f g (x,Pi;t), q = (u,d,s), q = (u,d,s), 
and Ji = J dpi/Ei. Note that we assign the label 1 to the momentum p and set 
f g (l)f g (2) / 9 (3)/ s (4)/ g (5) for loss and / s (4)/,(5) ^ / s (l)/ 9 (2)/ s (3) for gain in 
gg ggg process. We explicitly include the particle symmetry factor in the classical 
limit . 

We also obtain the evolution equations for quarks and antiquarks, 
p»dMx,p)= I [ I CW qq ^ qql [fMf q ,{A)-f q {l)f q ,{2)} 



+ 
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+ 



L L L \ [W ^ f^)fM) ~ W qq ^ 99 / g (l)/,(2)]. (3) 
/ / fc^[/,(3)/,(4)-/,(l)/,(2)] 

2 «/ 3 ^4 

i i ./• w *™v&W) - mm) 

L L L \ [W ^ f^)f g {A) - W qq ^ gg /,(l)/,(2)]. (4) 



where C — 1/2 if the final state consist of identical particles and C = 1 otherwise. We 
neglect in these Boltzmann equations those quantum mechanical effects coming from 
Bose enhancement factors (1 + f g ) for a gluon in the final state and Pauli blocking 
factors (1 — f q ) for final-state quarks or antiquarks. The quantum transition rates can 
be expressed in terms of differential cross sections, e.g. 



where s is the CM energy squared. 

These equations are highly non-linear and are impossible to solve analytically even 
with simplest initial state. To find solutions to these equations, we use the test particle 
method with Monte-Carlo sampling. 

3. Numerical Simulation 

The main idea to solve the Boltzmann equations (|^, |3], [|) is as follows. We assume that 
we know the initial state of the system, the position and momentum of each particle at 
a given time: 



We ignore in our study any collective phenomena which might be important, and we 
also ignore the quantum correlations and the nonperturbative properties of the strong 
interaction. With these simplifications, each particle moves freely as time proceeds, 
but collides occasionally with another particle if it comes close enough, i.e., within 



the distance determined by total cross section ^0^/71". In the collision, both particles 
abruptly change their momenta according to the differential cross section and/or produce 
new particles depending on the branching ratio of the various final states. The particles 
emerging from the collision are again moving freely until they make a collision with 
another particle. 

The Monte-Carlo simulation in general cascade code has several shortcomings: 

(i) The deterministic interpretation of the total cross section as a geometric criterion 
that decides whether two particles are colliding with each other or not is not a 
faithful representation of the quantum mechanical scattering process. It would be 




(5) 



A' 




(6) 
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more appropriate to use a collision probability as a function of the impact parameter 
based on a quantum transition amplitude. 

Another important issue is the assumption on the collision space-time: The two 
particles collide each other at their shortest distance or maximal force point (s). This 
assumption makes a big difference between a general (parton) cascade code and an 
simplified analytic solution scheme ||, |J or a relaxation time method |14) , 15 , [16 



The collision in a cascade code occurs only if there is substantial rapidity difference 
and/or transverse momentum difference since the shortest distance or maximal 
force point (s) between particles can be achieved in rather far future if two particles 
have same rapidity unless they have relatively high transverse momentum. On 
the other hand, most of collisions come from the small angle scattering between 
comoving partons in an analytic solution, namely collisions occurs among those in 
same rapidity. We believe this collision condition should be relaxed to allow those 
collisions between particles even if they are not in closest distance. 

We hope to address these issues in future publications. 
3.1. Initial Distribution 

It is very difficult to obtain a realistic initial distribution (^) for the simulation, because 
of the complexity of nuclear parton distributions and the dynamics governing the 
decoherence of the nuclear wavefunctions. The following choice is mainly governed 
by considerations of simplicity and easy implementation in the cascade code. Note that 
we are considering only head-on collisions, although our code is fully three-dimensional. 

As far as the transverse phase space distribution is concerned, reasonable choices 
can be derived from the semiclassical picture J7|. We assume that the number of 
produced partons is proportional to the number of primary collisions so that the 
transverse spatial distribution of the produced partons has the imprint of the transverse 
parton distribution of the colliding nuclei. This assumption can be expressed as the 
probability distribution for the initial transverse position (x,y) = (r± cos <p, r± sin</>): 

P(x,y) = B{l-^) (7) 

where R is the radius of the nuclei and B is a normalization constant. Once the 
transverse distance from the collision axis is selected, we can choose the azimuthal 
angle with equal probability over (0, 2ix). 

Krasnitz, Nara and Venugopalan (KNV) [17] found that the transverse momentum 
distribution from the study of small x-physics is given by the formula, 
1 dN 1 . / . \ 

/n(p±/A,) (8) 



where 



irR 2 dyd 2 p± g 2 



/ n (p±/A s ) = ai , p x < 3A S , 



exp( v /pl + m 2 /T eff) - i 
a 2 A 4 s H47ipJA s )pl\ P± >3A S (9) 
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where a± = 0.0295, a 2 = 0.0343, m = 0.067A S , T e ff = 0.93A S . y is rapidity and the 
saturation momentum A s = lGeV in RHIC and 2GeV in LHC. 

It is known that the rapidity of the produced particles after a heavy-ion collision is 
flat in the central rapidity region so that the rapidity y can be chosen from a constant 
distribution within the interval y m i n < y < y max - Once we have chosen the rapidity, the 
longitudinal momentum is given by the relation 



Pz — ym 2 + p^sinhy. (10) 

We further assume that all of the primary collisions occur in the reaction plane, 
z = 0, at time t = and the produced partons are born at the proper time, 
ro = \Jt 2 — z 2 , where tq is about 0.3 fm/c at RHIC and 0.13 fm/c at the LHC. We 
thus can find the initial time and longitudinal position of a produced parton once its 
rapidity is given, by calculating its velocity (3 Z = p z /E and setting f3 g = z/t to obtain 
t and z at the formation time ro. For simplicity, we initialize all partons at a common 
time t' = Tq and at z' = /3 z t . However, we do not allow the partons interact, until they 
reach a proper time r larger than the formation time ro. 

We assume that the number density per rapidity of produced partons is in the range 
dN/dy » 1000 ~ 1500. 

3.2. Cross Sections 

The dynamics enters into the Boltzmann equation through the total and differential cross 
section among partons. Differential cross sections at leading order a s for the processes 



are stated here explicitly for convenience, although they have been given elsewhere 18 



do-99^99 97^2 tu gu gt 



(3 - "5 " "ST - 75) (11) 



dt 2s 2 y s 2 t 2 u 

d a 99^qaq b na 2 .U t 9 t 2 + U 2 . . . 

—&r~ = 6* s *(t + u~ (12) 

da 9q ^ gq _ Ana 2 u s 9 s 2 + u 2 

~^^~^ 2 ~ { ~~s~u + l~^ 2 ~ ) (13) 



d<7«»»-9«» Ana 2 ,s 2 + u 2 . ,t 2 + s 2 2sl, 

— ST" = + 5 « b{ — ~ 3ut )] (14) 



4 fa 2 S 2 + M 2 t 2 + U 2 2U 2 
J t ~ -^[OacObd ~ 2 + OabOcd—^ abcd -—\ (15) 

do-0^99 32na 2 r u t 9t 2 + u 2 



dt ' 27s 2 QD 4 u 4 s 

and 

d a 9^99 da 9q ^ 9q da m ~* m da qq ^ qq 



Mt + T.- 7—7-] (I 6 ) 

(17) 



dt dt dt dt 
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We note that Eg. (|TT|, [12], fnj 0) are symmetric under the exchange of u «-> t as 
dictated by the indistinguishability of the identical particles. The classical cross sections 
must therefore include an overall factor 1/2. 

To get total cross sections, we integrate analytically the differential cross sections 
using the conditions, E cm > 1 GeV and the minimum momentum transfer p± > 0.5 
GeV so that the scattering angle 9 in the CM frame should satisfy the condition 
| sin^l > (0.5 GeV)/E, where E is the energy of the particle in CM frame. With this 
limitation on the scattering angle, we find c gg -, gg m 10/ GeV 2 and o~ qq ^ qq ~ 0.6/GeV 2 
and crg q ^g q ~ 0.7/GeV 2 for the strong coupling constant ot s = 0.3, which we will 
use throughout this work. These total cross sections are optimistic values for those 
perturbative processes because we choose low values for minimum CM energy as well 
as minimum momentum transfer. However, we note that about similar values for the 



cross section have been used by others [19 



We use the following approximation for the gg — > ggg cross section [10 
d a 99^999 _ 9C A a 3 s q\ 
dq\d y dk\~ 2 (q 2 ± + fi 2 D y 

Q(k±Xf — cosh?/)0( v /i — k± cosh?/) 



(18) 

where q± is the momentum transfer of two colliding gluons and k±_ the perpendicular 
momentum of the radiated gluon and y is its rapidity. Ca = 3 for the SU(3) gauge 
group, hd is a Debye screening mass, and A/ is the mean free path. The first step 
function represents the Landau- Pomeranchuk-Migdal(LPM) effect while the second one 
is the energy constraint of the radiated gluon. 

In order to integrate flT8|), we need several approximations, which we discuss now. 
For the Debye mass we use the value that is consistent with the gg — > gg cross section: 

A = 27rC ra a 2 /a™, (19) 

where C gg = 3. This approximation gives us hd ~ 0.4 GeV for o- gg ^ gg = 10/GeV 2 . The 
mean free path, A/, is defined as l/(p<Tr) where p is the density of particles and <jt is 
the average total cross section of the parton. This is a dynamical variable even if we 
fix the total cross section, because the system is rapidly expanding. We can estimate 
the initial density for heavy ion collision to be 0.34 ~ 0.5/GeV with 4000 ~ 6000 
gluons at r = 0.3 fm/c so that the mean free path is 0.2 ~ 0.3/GeV. But this number 
quickly increases as the system expands. As a reasonable approximation we choose 
A/ = 1.0/GeV in our study. With these two parameter choices, we can integrate the 
differential cross section numerically. Using the Runge-Kutta algorithm, we evaluate 
the total cross section a^ 9 ^ 999 ~ 3.6/GeV 2 

We can add the total cross section of each process to get the total cross section, 

e.g., 

a 9 ? = a 9 T 9 ^ 99 + a 9 ?-^ 999 + a 9 T 9 ^ ua + a 9 ^ + a 99 ^ sS . (20) 
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This total cross section will be used to decide whether two coming gluons will collide or 
not, by the criterion whether the impact parameter between two gluons in the reaction 
plane is less than ^Jorfn or not. Once it is determined that a collision takes place, we 
can use the branching ratio to select which channel will become active. After we choose 
the channel, we can use the differential cross section of the channel to select a scattering 
angle and momentum transfer. 

3.3. Basic Algorithm 

We first select the global laboratory frame as the reference frame in which the simulation 
is performed. We next consider all possible pairs of particles and calculate their impact 
parameter b. The impact parameter is defined as the shortest distance between two 
particles in their reaction frame, i.e. their mutual CM frame with one particle moving 
in the +z direction and the other moving in the — z direction. To determine b, we apply 
an appropriate Lorentz transformation to each pair of particles. Suppose the particle 1, 
which has the position (xi, yi, zi), moves along the +z-axis and particle 2, at (x 2 , 2/2, z 2)i 
moves along the — z-axis in the reaction frame. The impact parameter is explicitly given 
by: 

b = yj(x 1 -x 2 )* + (y 1 -y 2 )*. (21) 

We check whether a collision is possible by applying the criterion b < yar/TT as discussed 
before. We also check whether the condition z\ < Z2 is met. 

For a given particle, there maybe are many possible collisions but we keep only the 
earliest one, as observed in the laboratory frame. We next calculate where and when 
two particles will collide. In order to determine the collision point (s) for two particles, 
we make use of the "maximal force" principle, which is explained in detail in Appendix 
A. In brief, we postulate that the two particles are interacting through a retarded long- 
range field. This assumption is consistent with our use of lowest-order perturbative 
cross sections. We determine the point on each particle's trajectory where the retarded 
force becomes maximal. In general, these points correspond to different times in the 
laboratory frame, as an expression of the fact that the forces among particles propagate 
no faster than the speed of light. 

If the time difference between the two times corresponding to a single collision 
is greater than the mean free path, there is no collision between them, because they 
will be interrupted by neighboring particles. This can reduce the number of collisions 
substantially as the particle density goes up, or the relative CM energy of a particle pairs 
increases. This mechanism is similar to the procedure A of ref. [Ll] (restriction of the 
signal velocity). Here, however, we do not impose an ad-hoc cut on allowable collisions; 
rather, the finite signal velocity is embodied in the determination of the interaction 
points for each particle (see Appendix A). 

We then put all the collision events in time sequence and execute the collisions 
one by one in the global time frame. After each collision, the collision sequence must 
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be updated, because the particles emerging from collision have different phase space 
coordinates than before. 

3.4- Bulk Properties 

Since we have the full phase space information of all of particles composing our system, 
we can, in principle, calculate all the properties of the system. In particular, we can 
define the particle number current and the energy-momentum tensor as follows, 

N"(*) = f^jhrffrp), (22) 

T^(x) = J^p^f( x ,p). (23) 
The classical entropy current is defined as: 

= - j^f(x,p)\nf(x,p). (24) 

We also find it helpful to calculate the absolute momentum function weighted by the 
distribution 

P x (r,t) = [ d 3 p\p x \f(x,p), (25) 



P z (r,t) = j d 3 p\p z \f(x,p). (26) 

Local kinetic equilibrium requires P x = P z . 

The entropy of the system of classical particles cannot be calculated from (31) since 
the phase space of the test particles is given by a sum of 5-functions, 

f(x,p,t)= Y,S(x-*i(t))KP-Pi(t)), (27) 

i 

where is the position and momentum of the i-th particle. In order to define 

a smooth phase space density we smear this distribution in phase space such that 
the uncertainty relation is satisfied. We choose a spatial volume of extension Ax 
and a momentum space volume of extension Ap such that AxAp x > 1. We use the 
representations of the Dirac 5-function as a Gaussian: 

5(x - Xi) w /V a(x - X ' )2 (28) 

V 7T 

and smooth the probability up to a few nearest momentum bins and normalize the 
distribution so that the integrated particle probability remains unity. This is known as 
the coarse graining procedure and the entropy of a system is dependent on the procedure. 
In our calculation, we assign the probability, Eq.(^), up to the third nearest spatial 
bin from the position of the particle. Using this method, we calculate the one-particle 
entropy to be 3.2 for an isolated particle. 

4. Results and Discussions 



We are considering a system of which the total number of particles (gluons) is 6000. 
The particles have a flat rapidity distribution between +2 and —2. The total cross 
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Figure 1. Total number of partons as a function of time in 1/GcV unit: 6000 initial 
gluons and K — 2, 10, 20. Data are average over 20 runs which use different random 
number sets. 



section for the simulation has been calculated using the coupling constant a s = 0.3 and 
minimal momentum transfer 0.5 GeV as described before. We also include the effect of 
higher order diagrams in the form of a K- Factor K = 2. We use the KNV transverse 
momentum distribution (|]|D for the distribution. 

Figure 1 shows total number of particles as function of time. The primary partons 
are initialized at the proper time r (= 1.5/GeV at RHIC energy), and the secondary 
partons are produced by collisions between primary partons of the system. In the 
following, when we use the term "produced partons" , we refer to those partons that are 
created in collisions within our transport code, but not to the initialized partons that 
are produced before the start of our transport calculation. With our parameters, the 
secondary partons, mostly gluons, are not produced as abundantly as predicted in the 
analytical estimate of Baier, Mueller, Schiff and Son || . One reason for this difference 
is that the cross section (^) is small at early times due to the influence of the LPM 
effect, which provides for a strong suppression of radiation at high density. 

Figure 2 shows the ratio between transverse momentum and longitudinal 
momentum of the particles in a small box, which has the size of 2 x 2 x 1 fm 3 and 
is located at the coordinate origin, as a function of time. This ratio should be close 
to unity when the kinetic equilibration is reached. The figure shows that the ratio is 
increasing beyond unity, implying that the system is thinning more rapidly along the 
collision axis than in the transverse direction for K = 2 and even for K = 10. Kinetic 
equilibration can be achieved for the unrealistically large cross section with K = 20. 

Figure 3 shows the energy density of the box as a function of time. The high energy 
particles escape from the box very quickly and a relatively slow expansion follows. 

Figure 4 shows the particle number density of the box which we are considering 
as a function of time. We can see that more and more particles are staying in the box 
as the cross section gets larger, but the densities are almost the same after those high 
momentum particles escape. Scattering occurs only occasionally as time goes on. 

In conclusion, we do not observe kinetic equilibration for a realistic cross section 
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Figure 2. |Px|/X] lf>zl as a function of time in 1/GeV unit: 6000 initial gluons and 
K-Factor = 2(+), 10(x) and 20(*). Data are the average over 20 runs in which each 
run uses different set of random numbers. A small box (2 x 2 x 1/m 3 ) at center is used 
for measurement. 
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Figure 3. The energy density per fm 3 in the box at central region as a function of 
time . 
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Figure 4. The number density per fm 3 in the box at central region as a function of 
time . 
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in our simulations. However, when we increase the cross section substantially (10 — 20 
times!) approximate isotropy of the momentum distribution can be achieved. In our 
investigation here, we kept the low-momentum fixed. On the other hand, dynamical 
screening in an expanding medium provides an effetive low-momentum cutoff that varies 
with time. This raises the question whether a selfconsistent treatment, in which the 
time-dependent density is used to determine the time variation of the total cross section, 
can exhibit an approach to equilibrium without unrealistic assumptions. We hope to 
return to this question, as well as to the issue of including 3 — ^ 2 processes to ensure 
detailed balance, in a future publication. 
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Appendix A. Collision points 

Quantum mechanics does not allow us to determine the precise moment of when 
a collision between two particles occurs due to the uncertainty relation. Here, we 
discuss this question in the context of classical mechanics, in particular, in classical 
electrodynamics. Consider two particles with charges q±, q2, and positions and velocities, 
m m ) and (r",v v ), respectively. The electromagnetic field induced by particle 2 at 
the position of particle 1 is given by ]20 



F ^ = [v {x -r { r' I W [(X - r(<))V (A1) 
-(i-rte))V] 

where r{r' x ) is the source point associated with the field at point x. Each particle moves 
along its trajectory at constant velocity: 

x" = + u^t, (A.2) 

r» = W + uV. (A.3) 

Using the retardation condition, we can obtain the proper time of the source point for 
the field point x: 

r> x = ( x -b)-v- ^[(x-b)-v] 2 -(x-b) 2 . (A.4) 

The force on particle 1 due to the field of particle 2 is given by 

e[(x - r'Yv ■ u - (x - r') ■ uf] 

v [[(a - b + ur) ■ v] 2 -{a-b + ur) 2 fl 2 1 ' 
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where r' = r(r' x ). This force will be maximal when the denominator is minimal. This 

condition will be met at 

(a - b) ■ v(u ■ v) - (a - b) ■ u 
t = 2 7 \2 • ( A - b > 

Even though the classical collision occurs gradually, we associate this space-time point 
with the collision moment of the particle 1: 

x" = a" + «"r . (A.7) 

Using the same procedure, we obtain the maximal force on particle 2 by the particle 1 
at the proper time 

/ (6 - a) • u(v ■u)-(b-a)-v 
r = (A. 8) 

v z — (u • vy 

and the collision space-time point 

r" = 6" + u"t£. (A.9) 

Since the definitions of r and Tq are manifestly Lorentz invariant, the interaction points 
are independent of the reference frame in which the two-body collision is calculated. 

The force between the two particle is more complex due to the presence of color 
degrees of freedom. However, in lowest-order perturbation theory, the space-time 
behavior of the QCD force is identical to that in electrodynamics. We can thus apply 
the equations derived above. 

Appendix B. Monte Carlo Sampling 

We explain here the Monte-Carlo sampling for the process gg — > gg only. The other 

processes are the same except the process gg — > ggg which we discuss briefly at the 

end. We change the variable from the Mandelstam t to rj — cos 9. The differential cross 

section becomes 

da^gg _ Q na 2 (i _ 7? )(i + ^) 

where 77 is confined by the minimal momentum transfer Q which we will choose either 
0.5 or lGeV. The total cross section can be obtained analytically by integrating over 
the allowed range. The natural way to choose 77 by monte Carlo method is that we 
calculate 

r= / — (B.2) 

and invert this equation to solve for 77. Since this is not easy to do, we divide the 
(Vmin,Vmax) m to several intervals. Each interval can be integrated over analytically to 
give an area. We generate a random number 7*1 and select the interval depending on the 
area, i.e., if r x < si/a T , the interval 1 is chosen, and so forth. Once we have selected 
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the interval (r] start ,r] end ), we can find the straight line connecting the points at rj start 
and r} en( i and choose r\ between (r] start , rj end ) by linear random sampling and then use the 
accept-reject method. This triple Monte Carlo sampling gives a satisfactory result. 

The sampling of the gg — ^ ggg process is slightly more complicated. Using the step 
function in (|I8D, we find the rapidity of the bremsstrahlung gluon, (y m in, Umax), where 
for a given CM energy y/s, 



y max = In \J y/sXf + \J y/sXf - I J . (B.3) 
Once we have chosen the rapidity, we have 

(B.4) 




and 

< q\ < s. (B.5) 

Again, we can use the accept-reject method to sample (k"j_, qj_). This algorithm generates 
an acceptable distribution within a reasonable CPU time. 
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